classdef GetMagnetometerSensitivity < Table1.NMRGJobs.NMRGJob
    methods
        function obj = GetMagnetometerSensitivity(nmrg)
            obj@Table1.NMRGJobs.NMRGJob(nmrg, 'GetMagnetometerSensitivity');
            obj.defineParameter();
        end
        
        function obj = defineParameter(obj)
            obj.paramParser=inputParser;
            obj.paramParser.addParameter('nCol',             1024, @(x) x>0);
            obj.paramParser.addParameter('average',          100, @(x) x>0);
            obj.paramParser.addParameter('signal_frequency', 30.0, @(x) x>0 && x < 100);
            obj.paramParser.addParameter('signal_amplitude', 0.02, @(x) x>0 && x < 0.1);
            obj.paramParser.addParameter('isReport',         false, @islogical);
            obj.paramParser.addParameter('description',      'No description.', @ischar);
        end
        
        %%
        function obj = preprocess(obj)
            param = obj.dataPackage.parameters;
            nmrg = obj.boss;

            nmrg.setParametricMagnetometer('frequencyZ', nmrg.SystemInfo.WorkingPoint.mainFrequency, ...
                                           'isLockTransverseField', true);
            nmrg.bx.setAC('frequency', param.signal_frequency, 'amplitude', param.signal_amplitude);
        end
        
        %%
        function obj = dataAcquisition(obj)
            import MeasurementComponent.Signal.*

            nmrg = obj.boss;
            param = obj.dataPackage.parameters;
            measurement_data.start = datetime('now');
            
            nmrg.bTrans.unsubscribe_signal();
            nmrg.bTrans.subscribe_signal(ziDemodSpectrum.X_fft_abs_avg);
            measurement_data.spectrum = nmrg.bTrans.getDemodSpectrum(param.nCol, ...
                                                                    'average', param.average, ...
                                                                    'fSpan', nmrg.bTrans.getSetting('bandwidth'));
            

            measurement_data.finish = datetime('now');
            obj.dataPackage.measurement_data = measurement_data;
        end
        
        %%
        function obj = analyze(obj)
            measuredData = obj.dataPackage.measurement_data;
            signal = measuredData.spectrum;
            nmrg = obj.boss;

            %------------------------
            V2nT = nmrg.bx.setting.volt2curr * nmrg.bx.setting.curr2nT; % control voltage [V] to actual field [nT]
            slope = nmrg.SystemInfo.WorkingPoint.slopeX;                % control voltage [V] to signal amplitude [V]
            coeff = abs(slope/V2nT);                                    % [V/nT] signal amplitude generated by 1nT transverse field in x direction

            enbw = signal.getENBW(nmrg.bTrans.subscribed_signal{1, [1 2]});
            data = signal.getNodeData(1);
            psd=data.y/sqrt(enbw) / coeff * 1e3;
            [ncount, edges] = histcounts(psd, 0:0.1:100); % count sensitivity fro 0 ~ 1000 pT/sqrt(Hz)
            [~, loc, wid ] = findpeaks(ncount, 'NPeaks', 1, 'SortStr', 'descend');

            %------------------------
            analysis_result.enbw  = enbw;
            analysis_result.V2nT  = V2nT;
            analysis_result.slope = slope;
            analysis_result.coeff = coeff;
            analysis_result.meanSensitivity = mean( psd( psd<edges(ceil(loc+3*wid)) ) );
            analysis_result.fList = data.x;
            analysis_result.psd   = psd;
                        
            obj.dataPackage.analysis_result = analysis_result;
        end
        
        %%
        function obj = postprocess(obj)
            param = obj.dataPackage.parameters;
            measuredData = obj.dataPackage.measurement_data;
            analysis = obj.dataPackage.analysis_result;
            nmrg = obj.boss;
            
            % Post operations
            nmrg.bx.disable();
            
            
            % summary
            summary.start = datestr(measuredData.start);
            summary.amplitude = param.signal_amplitude * analysis.V2nT;
            summary.frequency = param.signal_frequency;
            summary.psdBackground = analysis.meanSensitivity;            
            summary.finish = datestr(measuredData.finish);
            
            obj.dataPackage.summary = summary;
            nmrg.SystemInfo.MagnetometorSensitivity = summary;
            
            % Store Data
            getSession;
            sess.addData('GetMagnetometerSensitivity', obj.dataPackage, param.description);
        end

        %%
        function obj = presentation(obj, viewer)
            param        = obj.dataPackage.parameters;
            measuredData = obj.dataPackage.measurement_data;
            analysis     = obj.dataPackage.analysis_result;
            summary      = obj.dataPackage.summary;
            
            if nargin == 1
                f = figure('Name', sprintf('Sensitivity[%s]',datestr(measuredData.start)), 'WindowStyle', 'docked');
                ax = subplot(1, 1, 1, 'Parent', f);
            else
                Table1.NMRGJobs.NMRGJob.updateViewer(viewer, param, summary);
                ax = subplot(1, 1, 1, 'Parent', viewer.gui.uiHandle.panelPlot);
            end
            
            
            fList = analysis.fList;
            psd   = analysis.psd;
            enbw  = analysis.enbw;
            sig_amp  = summary.amplitude;
            sig_freq = summary.frequency;
            meanSensitivity = analysis.meanSensitivity;
            
            semilogy(ax, fList, psd, 'b.-', [fList(1), fList(end)], [meanSensitivity, meanSensitivity], 'r--');
            xlabel(ax, 'frequency (Hz)'); ylabel(ax, 'sensitivity pT/sqrt(Hz)');
            xlim(ax, minmax(fList(:)'));
            title(ax, sprintf('Signal = %3.2f nT @ %2.1f Hz, ENBW = %3.2f Hz', sig_amp, sig_freq, enbw));
            legend(ax, {'PSD', sprintf('Background = %3.2f pT/\\surd{Hz}', meanSensitivity)}, 'Location', 'best', 'Interpreter', 'tex');
            grid(ax, 'on');
        end
    end
    %%
    methods(Static)
        function res = compare(objList, varargin)
            p = inputParser;
            p.addParameter('isPlot', true, @islogical);
            p.addParameter('viewer', [], @(x)  isa(x, 'Table1.GUI.DataViewer') || isempty(x));
            p.parse(varargin{:});
            
            res = compare@Table1.NMRGJobs.NMRGJob(objList, 'viewer', p.Results.viewer);

            if p.Results.isPlot
                startTime = arrayfun(@(t) t.start, res.summary, 'UniformOutput', false);
                if isempty(p.Results.viewer)
                    f = figure('Name', 'Sensitivity');
                    ax = subplot(1, 1, 1, 'Parent', f);
                else
                    ax = subplot(1, 1, 1, 'Parent', p.Results.viewer.gui.uiHandle.panelPlot);
                end
                fList = cell2mat(arrayfun(@(d) d.fList', res.analysis, 'UniformOutput', false));
                psd = cell2mat(arrayfun(@(d) d.psd', res.analysis, 'UniformOutput', false));
                
                semilogy(ax, fList, psd, '.-');
                xlabel(ax, 'frequency (Hz)'); ylabel(ax, 'sensitivity pT/sqrt(Hz)');
                xlim(ax, minmax(fList(:)'));
                legend(ax, startTime, 'Location', 'bestoutside', 'Interpreter', 'tex');
                grid(ax, 'on');
            end
        end
    end
end

